Transmission dynamics and baseline epidemiological parameter estimates of Coronavirus disease 2019 pre-vaccination: Davao City, Philippines

The Coronavirus disease 2019 (COVID-19) has exposed many systemic vulnerabilities in many countries’ health system, disaster preparedness, and adequate response capabilities. With the early lack of data and information about the virus and the many differing local-specific factors contributing to its transmission, managing its spread had been challenging. The current work presents a modified Susceptible-Exposed-Infectious-Recovered compartmental model incorporating intervention protocols during different community quarantine periods. The COVID-19 reported cases before the vaccine rollout in Davao City, Philippines, are utilized to obtain baseline values for key epidemiologic model parameters. The probable secondary infections (i.e., time-varying reproduction number) among other epidemiological indicators were computed. Results show that the cases in Davao City were driven by the transmission rates, positivity proportion, latency period, and the number of severely symptomatic patients. This paper provides qualitative insights into the transmission dynamics of COVID-19 along with the government’s implemented intervention protocols. Furthermore, this modeling framework could be used for decision support, policy making, and system development for the current and future pandemics.


Introduction
The Coronavirus disease 2019  caused by SARS-CoV-2 virus has globally spread since its first detection in Wuhan, China, last December 2019 and has caused millions of deaths worldwide [1][2][3][4]. From then on, different control measures have been implemented to mitigate its spread. In the absence of vaccines, different countries have resulted in the implementation of non-pharmaceutical interventions (NPIs) [5][6][7][8][9][10], such as the quarantine measures/isolation, contact tracing, physical distancing, and mass testing. In the Philippines, various levels of community quarantines have been implemented in different parts of the country since the first confirmed COVID-19 case was reported [11]. These community quarantines are classified as enhanced community quarantine (ECQ), modified enhanced community quarantine (MECQ), general community quarantine (GCQ), or modified general community quarantine (MGCQ). Regardless of the quarantine classification, the general public is urged to stay in their respective households, practice physical distancing, and wear face masks and face shields to lessen the spread of the disease. School closures, halting of mass gatherings and other potentially superspreading events (e.g., religious activities) have been imposed. NPIs are continually being implemented to prevent a surge of cases due to the crowding of individuals [12,13]. Noteworthy, however, is the apparent difference in the strictness in implementing different control measures between quarantine classifications as summarized in Table 1 [14]. Furthermore, the number of deployed quarantine enforcers and the number of monitoring stations put up differs among respective quarantine classifications. Even between regions in the Philippines under similar quarantine classification, differences in resources and workforce, among others, contribute to variability in the efficacy of the intervention. Hence, assessing the effects of the different quarantine classifications on a regional, provincial, or city-level could help decision-makers formulate a more appropriate policy to mitigate the health crisis.
A baseline measure should be developed to properly comprehend the effect of an intervention to the COVID-19 cases; this will be a basis for determining any subsequent changes brought about by these interventions [15,16]. Epidemiological models help identify critical factors of non-pharmaceutical interventions (NPIs) implementation, such as timing and frequency, that support the control of disease spread [17]. Understanding the parameters that determine the course of an epidemic is critical for health-related decision-making, as it allows for the development of disease mitigation and control methods, as well as the provision of treatment to individuals who have been infected or become ill [18]. These models have been used as guides to policy and decision-makers as well as implementers to combat outbreaks [19][20][21]. To this end, several mathematical models on the dynamics of COVID-19 with NPIs have been developed and published to project the different COVID-19 transmission scenarios  [22,23]. Nevertheless, assessing the effects of different quarantine classifications to key epidemiological parameters on a more specific scale, i.e., regional level, has not been widely documented. The high transmissibility and virulence of SARS-CoV-2 resulted in a significant number of severe and critical cases requiring specialized treatment and intensive care beds-forcing the development of predictive models capable of estimating healthcare demands and assisting decision-making [24][25][26]. For better contextualization of the model, local COVID-19 dynamics and NPIs in Davao City were used in the current case study. Davao City is the largest city in the Philippines based on land area. It is located in Region 11 (Davao Region), the most populous region in Mindanao-the southern part of the country [27]. The statistics from the Department of Health-Davao Center for Health Development (DOH-DCHD) show that Davao City accounts for almost 57% of Davao Region's COVID-19 cases (as of June 26, 2021). Davao City is also the capital of Region 11, serving as the focal business hub and other activities in the region. Hence, the importation and transmission risk of COVID-19 from and to neighboring provinces and regions is high.
This paper modified the classic Susceptible-Exposed-Infected-Recovered (SEIR) deterministic compartmental model [28] into a Susceptible-Exposed-Infectious Hospitalized-Infectious Monitored-Recovered (SEI H I M R) model to describe the dynamics of COVID-19 in Davao City incorporating the Philippines' Department of Health's quarantine and isolation protocols, namely quarantining exposed individuals and isolating confirmed positive cases according to the severity of symptoms or presence of comorbidity [13]. This model incorporated an outflow from the E compartment back to S similar to [29]-a dynamic that has not been widely applied to models in a Philippine regional setting. This link signified that, instead of the standard classification of the E compartment [28], not all exposed individuals in this model are on a latent stage of the disease, do not become infectious, and hence become susceptible again after they were released from quarantine or isolation. Moreover, the model also incorporates the progression or worsening of symptoms of some previously asymptomatic positive cases (i.e., presymptomatic). Disease-related deaths also in this model are assumed to only occur among severely symptomatic patients. Furthermore, in Davao City's context, isolation and, or quarantining of the infectious individuals do not guarantee that the patient has stopped contributing to the pathogen's transmission; viral spread may occur before an intervention due to delays in test result. Moreover, health and emergency workers can still get infection from quarantined or isolated individuals. Hence in this model, we assumed that the compartments I H and I M equally drive the transmission regardless of isolation or quarantine. Nevertheless, the contact tracing method in Davao City during the period covered in this paper is assumed to be efficient whereby minimizing the possibility of an untraced infectious individual. Thus, the model does not incorporate a compartment for undetected cases to make analysis tractable.
We aimed to estimate the baseline epidemiological parameters, i.e., pre-vaccination COVID-19 period, using the least-square method and bootstrapping techniques in quantifying parameter uncertainty. Epidemiological measures such as the basic reproduction number(R 0 ), statistical time-varying reproduction number (R s t ), and the deterministic time-varying reproduction number (R d t ) have also been computed. These numbers generate insightful thresholds on the secondary infections generated by a COVID-19 infectious patient upon interaction with the susceptible population [30]. We also provided basic local stability analysis of the model. The COVID-19 outbreak has led researchers to investigate numerous factors of disease transmission and evolution [31,32]. The empirical results from the generated model can guide the local government unit in reviewing their respective protocols to mitigate the spread of COVID-19. Furthermore, the model described here is general enough to be used in studying COVID-19 dynamics in other regions, cities, or municipalities.

Ethics statement
In compliance with the Joint Memorandum Circular No. 2020-0002 by the Philippines' Department of Health on the data privacy guideline, processing and disclosure of COVID-19-related data for disease surveillance and response, the study was reviewed by the Department of Health XI Cluster Ethics Review Committee with protocol number P211111601. Furthermore, this study does not involve human participants or personally identifiable information.

Model formulation
In this paper, the closed homogeneous population was divided into five population classes, namely, susceptible (S), quarantined/exposed (E), infectious individuals treated in either a monitoring facility (I M ) or hospitals handling severe cases (I H ), and recovered (R) individuals. We specifically defined the compartments as follows: • Susceptible (S): Individuals who have not been infected with COVID-19 but are at risk of contracting the disease; • Exposed (E): Individuals who were under quarantine who may have been exposed to the disease and were still in the latent stage or who may have been exposed but have not necessarily contracted the virus;  Table 2. Herein, the mean residence time in E, in I M , and in I H are respectively given as The system of differential equations that represents the transmission dynamics of the SARS-CoV-2 virus in the city was obtained by adding the inflow rate and subtracting the outflow rate for each respective component. The following system (1) mathematically described the disease dynamics

with/where S(t), E(t), I M (t), I H (t), and R(t) are nonegative, and N(t) = S(t) + E(t) + I M (t) + I H (t) + R(t).
The parameters were piecewise functions of time subject to the duration of implementation of the different levels of quarantine classifications (i.e., ECQ, MECQ, GCQ, and MGCQ). It could be easily seen from the system (1) that because Hence, the model is well posed, and all solutions remain in O.

Data
Publicly available data on the daily COVID-19 new infected cases in the Philippines can be accessed through the government data drop [35]. However, the data used in this study were obtained under a partnership and non-disclosure agreement from the DOH-DCHD: Regional Epidemiology and Surveillance Unit (RESU). This partnership and data-sharing agreement were sought for better and more efficient data, and result validation, which would have been difficult if the data was acquired from the data drop. The epidemiological dataset used included dates of illness' onset, dates of surveillance report, and the health status upon admission of confirmed cases (e.g., asymptomatic, mild, moderate, severe, and critical). The dataset was between March 8, 2020, to March 5, 2021. We categorized the epidemic data according to two groups: the infected individuals under monitoring in TTMFs (e.g., asymptomatic and mild cases), and the infected patients admitted to COVID-19 referral hospitals (e.g., moderate, severe, and critical cases). Even though the working model is a baseline model that does not distinguish vaccination as a disease control measure, the model could be used to assess pre-vaccination transmission dynamics. Hence, we only used data before the start of the vaccine rollout in Davao City, which was on March 5, 2021. For reproducibility, a sample dataset which is in compliance with the non-disclosure agreement could be accessed through [36].

Community quarantine timeline pre-vaccination
Pre-vaccination, the city government of Davao, following the recommendations and guidelines issued by the Inter-Agency Task Force (IATF), imposed NPIs such as social distancing, lockdowns, curfews, and office and school closures, among others, to reduce disease transmission. The intensity of the implementation of NPIs varied according to the type of community quarantine imposed as summarized in Table 1. Thus, we divided the epidemiological data into four distinct periods according to the intensity and classification of community quarantine implemented in the city.

Sensitivity analysis
Sensitivity analysis was performed a priori to parameter estimation. This was used as the basis why certain parameters needed to be identified while some were fixed. Parameters with higher sensitivity needed to be estimated reliably subject to the available information. We performed uncertainty and sensitivity analysis using the Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficient (PRCC) method [37]. To perform the LHS analysis, 11 of the 13 system parameters plus a dummy were varied simultaneously while we fixed the values of the α and δ parameters according to

Model parameterization
Model parameters were estimated following the methods of Chowell, and Banks et al. [39,40]. The model was calibrated to the daily new infected cases in Davao City from March 8, 2020 (the first case of COVID-19) to March 5, 2021. The model calibration was performed individually for each quarantine period considered (e.g., Q1-Q4). The model has 13 epidemiological parameters. Two (2) of these parameters were estimated: the transmission rate (β), and the positivity proportion (r). First, a seven (7)-day moving average filter was applied to smoothen the random variations in the daily incidence data. To fit the asymptomatic-or-mild daily incidence data to the model, we defined the newly infected people under monitoring as The moderate-to-critical daily incidence data were fitted to the model using the new hospitalized infected people defined as H(t) = δ H rqE. The data-fitting problem was then solved using the least squares (LS) technique given by the objective function f^t i Þ ¼ argminS n i¼1 ðf ðt i Þ À yðt i ÞÞ 2 , y(t i ) was the observed daily new cases, and f(t i ) was the corresponding model simulation.
In particular, f(t i ) = (M(t i ), H(t i )) and y(t i )'s are the daily reported monitored and hospitalized cases. The fitting procedure was executed using the MATLAB lscurvefit function with numerical optimization through a trust-region reflective algorithm [41]. Other parameter values were either estimated from the data [42,43] or taken from existing published literature. The complete list of the system parameters and their sources is presented in Table 2. The estimated parameter values are presented in Table 3. These values were then analyzed for feasibility and were verified or compared to other existing literature. The analysis is expounded in the results section.

PLOS ONE
Transmission dynamics, baseline epidemiological estimates of COVID-19 pre-vaccination

Bootstrapping method
The bootstrapping method was used to simulate the lower and upper bounds (95% confidence level) around the model fit in assessing parameter identifiability [39]. The parametric bootstrapping approach generated 10,000 bootstrap realizations assuming a Poisson error structure of the best fit from the LS fitting method. The generated bootstrap realizations were used to resimulate the SEI H I M R model and to derive the empirical distributions of the estimated parameters within a 95% confidence interval.

Time-varying reproduction number
The daily statistical time-varying reproduction number (R s t ) in this paper was computed following the methods of Cori and colleagues with prior mean μ SI = 2.6 and standard deviation σ SI = 2 [37,44] and mean μ SI = 4.8 and standard deviation σ SI = 2.3 [45]. Furthermore, the COVID-19 cases according to the onset of symptoms, were used to statistically estimate the R s t . On the other hand, a formula for the deterministic time-varying reproduction number (R d t ) was also derived from the compartmental model using the next-generation matrix method [30]. The R s t was then compared to the R d t .

Qualitative analyses
Some of the qualitative analyses performed in this study involved computing the basic (R 0 ) and time-varying deterministic reproduction numbers (R d t ), and solving for the model's local stability, equilibrium points, and elasticity indices. The detailed steps to these analyses, proof and calculation can be found in the supplementary S1 File. The basic reproduction number (R 0 ) was the initial reproduction number of the pathogen at the start of the pandemic, whereas the deterministic time-varying reproduction number R d t was the reproduction number at any particular point in time t. The R d t described the average number of secondary infections caused by an infective individual in a susceptible population in the presence of control measures [38]. It provided a threshold for disease outbreak: if R d t < 1, the disease cannot persist; when R d t � 1, the virus can spread within the population where each infectious individual produced at least one new infection [46,47]. The derived closed-form formulas of R 0 and R d t were as follows: To further determine how the parameters affect R 0 , we solved for the elasticity index (normalized sensitivity index) [48]. The elasticity index measures the relative change of R 0 with respect to a parameter, for example, β, denoted by U The following are the proportionality and the elasticity indices for the relevant parameters: 1. The direct proportionality of the transmission rate β to R 0 is supported by its elasticity index of 1. This means that a unit increase of β results in a unit increase of R 0 . Hence, any increase or decrease of the transmission of the pathogen from one person to another directly affects the size of the pandemic and its mitigation. Conversely, this result supports that the measures designed to hamper this transmission, such as wearing facemasks, hand washing, disinfection, observing physical distancing, and isolation, among others, [49,50] can directly slow down the spread of the disease.
2. The direct proportionality of the positivity proportion r to R 0 is also quantified with an elasticity index of 1. This means that a unit increase in r also results in a unit increase of R 0 . Thus, interventions designed to decrease r could directly lead to the mitigation of the spread of the disease. The interventions related to r are, but not limited to, contact tracing and testing. An increase in the effort in these interventions will directly increase the identified positive cases, thereby helps curb the pandemic for better management. That is, the more contact traced, quarantined, and tested individuals, the better, regardless if only a few among them become a confirmed positive case. It is even more ideal to have a lower positivity rate among tested individuals.
3. The elasticity index of A on R 0 is −1. This backs the inverse proportionality of A to R 0 , which means that a unit increase in A results to a unit decrease of R 0 . Conversely, as A is just the inverse of the mean residence time in E (A −1 ), a unit increase/decrease in A −1 would contribute to a unit increase/decrease in R 0 , i.e., the faster the individual transitions out of the E compartment (smaller residence time), the more manageable the spread of the pathogen will be (with a reduced R 0 ). Since each term of A is an outflow rate identified with an intervention, this result lends support to more immediate COVID-19 test results combined with a reasonably shorter quarantine/isolation period (days) for reducing the risk of viral spread. However, due to the biological dynamics of COVID-19, there is a minimum required quarantine/isolation period, which imposes a positive lower bound for practical values of A when isolation is combined with limited testing capacity.
4. The inverse proportionality of B to R 0 is corroborated by its elasticity index of −0.99937. A unit increase of B would lead to a decrease in R 0 by 0.99937 units. Conversely, a unit increase in the mean residence time in the I M compartment (B −1 ) would lead to an increase of 0.99937 units to R 0 . This signifies that the more number of asymptomatic and, or mildly symptomatic COVID-19 cases, the more the pandemic can be controlled. Hence, measures in ensuring that the vulnerable population gets due protection and medical attention produce positive results in preventing the further spread of the virus and disease-related deaths [51,52].
The disease free, and endemic equilibrium points of the model were also solved and were verified by the proof of its local stability. This meant that two scenarios could be expected to come about in the COVID-19 pandemic dependent on certain conditions: (1) the disease would die out, or (2) becomes endemic [48]. These conditions were then verified to Davao City's circumstances, for which it was found that endemicity would, at best, be the more favorable scenario to be expected, as the disease free scenario is relatively too difficult to achieve, and at the point when this paper is written, is not yet realistic. Hence, policies toward the management of COVID-19 should consider an endemic scenario where the disease continues to exist but on a controlled manner. Recommendations relative to this are further expounded in the Conclusion part of this paper.

Sensitivity analysis
The sensitivity analysis helped ascertain which parameters were to be identified or which parameters were to be fixed. Fig 2 presents the LHS/PRCC analysis results for the different response functions considered with the corresponding values presented in Table 4. It was observed that the parameters β and r were strongly positively correlated to I H and I M cases. In addition to β and r parameters, δ M was also strongly positively correlated to I M cases. These indicated a strong association between high values of the β, r, δ M parameters to the number of COVID-19 cases, specifically during the early stages of the pandemic. The correlation of these parameters to the response functions decreased over time during the course of the pandemic. This result supported the need to distinguish these three (3) parameters as essential factors in the transmission dynamics of the virus. However, as δ M is dependent on the biology of SARS-CoV-2, intervention recommendations focused mainly on policy-preventable parameters β and r. Moreover, as the severity proportion q and disease-induced death rate μ H directly equates to COVID-19 morbidities that we also put attention to these parameters. The implementation of control interventions were essential for the reduction of the values of these parameters (β, r, q, and μ H ) to mitigate the spread and lessen the hospitalizations due to COVID-19. Logrosa and colleagues corroborate that if the local government prioritized reducing COVID-19 fatality, the COVID-19 pandemic will be manageable, at least in the context of Davao City, Philippines [27].

PLOS ONE
Transmission dynamics, baseline epidemiological estimates of COVID-19 pre-vaccination

Epidemiological parameters
From the gathered data, Fig 3 shows the COVID-19 daily and cumulative incidence in Davao City over the date of illness onset and the four-time periods that cover the community quarantine levels from March 8, 2020, to March 5, 2021. The complete list of system parameters as a result of the data fitting is presented in Table 3. Fig 3 showed the observed COVID-19 cases and the best model fit to the data and that the model was able to approximate the dynamics of

PLOS ONE
Transmission dynamics, baseline epidemiological estimates of COVID-19 pre-vaccination the pathogen. Following the bootstrapping method, the resulting empirical distribution of parameters β and r, including R d t and the 95% confidence intervals of the estimated parameter values are presented in Fig 4. These show that the distribution is normal and that the model and the best fit values were able to capture the transmission dynamics of COVID-19 in Davao City, Philippines.
It was observed in the Figs 3 and 4 that the estimated values of the transmission rate (β) were at its highest at the start of the pandemic and increased as the quarantine lockdowns were eased over time (compare Table 3 [53], where the strictest measure implemented in ECQ (Q1) led to the least transmission rate at the next classification (Q2). However, since right after Q1 (ECQ; strictest), Q2 (GCQ; most permissive) was implemented, its perceived lag effect brought about a higher transmission rate in Q3, which was almost twice as much as in Q2. The abrupt relaxation of quarantine classification from the strictest to the most relaxed one led to high transmissibility of the disease as more people were allowed mobility. Hence, it would have been best if a gradual relaxation of protocols were followed instead. Nevertheless, comparatively, these estimated values are similar and or within the range of published estimates throughout the world [22], for example, in India [51], and China [52]. Hence, the same as the case of China, with proper stringent preventive measures, the COVID-19 cases in Davao can still be brought down to manageable levels. Table 3 show that the model estimates of the positivity proportion r were lowest during Q1 and highest during Q2. The results mean that 89.87% of the identified "close contacts" did not develop the disease during Q1, 82.42% during Q2, 82.98% during Q3, and 87.62% during Q4. As Davao City had strict mobility restrictions during Q1 may have led to lower contact among individuals and may have lessen the transmission of the pathogen, thus the tracing and testing capacities were more efficient hence the lowest r value. However, as lockdowns were eased in Q2 and Q3, we observed increases in r despite the tracing and testing efforts. This is due to the

PLOS ONE
Transmission dynamics, baseline epidemiological estimates of COVID-19 pre-vaccination increase of the interacting population size while the tracing and testing capacity remained constant. Furthermore, the testing and tracing guidelines during Q2 and Q3 were a bit more relaxed causing a fewer identified "close contacts", hence the increase of r. The decrease of these values would have been a favorable outcome towards controlling the spread of the disease [27].
Different scenarios were considered concerning the different r values estimated as shown in Fig 5. It shows that if mass testing was done to at least halve the positivity rate, the daily new exposed population would have reduced from 73 to 21 individuals after 60 days of Q1 and from 34 to 14 individuals during Q2; the daily new I M cases would have reduced from 24 to 5 individuals after 60 days of Q1 and from 30 to 8 individuals during Q2; and the daily new I H cases would have become 0 from 2 individuals after 60 days of Q1 and 0 indv from 1 individuals during Q2. Moreover, if the efforts were a bit laxer to twice than what was the positivity rate, the daily new exposed population would have increased to 358 individuals after 60 days of Q1 and to 90 individuals during Q2, the daily new I M cases would have increased to 175 individuals after 60 days of Q1 and to 122 individuals during Q2, while the daily new I H cases would have become 13 individuals after 60 days of Q1 and 3 individuals during Q2. These results corroborate with the need for mass testing, contact tracing, and the need for faster release of test results towards mitigating the disease [54].
With respect to the severity proportion q as shown in Table 3, we found that 77.29% (1 −q = 0.7729) of the confirmed cases are asymptomatic to mild, while 22.71% (q = 0.2271) were moderate to critical cases during Q1. Likewise, during Q2, 92.12% are asymptomatic to mild cases, and 7.88% are moderate to critical cases. The moderate to critical cases during Q3 slightly increased with 8.42% of the confirmed cases and 7.42% moderate to critical cases during Q4. The majority of the confirmed cases in Davao City consisted of asymptomatic to mild cases. However, even though the ratio of moderate to critical cases in Davao to asymptomatic to mild cases is relatively low, its total number of cases still amounts to more than Davao City's

PLOS ONE
Transmission dynamics, baseline epidemiological estimates of COVID-19 pre-vaccination total COVID-19 bed capacity. Hence, Davao COVID-19 referral hospitals are usually maxed out of their capacities [55]. Policies and budgetary support towards the increase of Davao City's health system and increasing hospital capacities should be prioritized.  [55].
As seen in Fig 6, the overall trend of the deterministic reproduction number R d t follows that of the statistical reproduction number R s t . Furthermore, there is a decrease in the time-varying reproduction number during the Q2 period as compared to Q1. This is during the early stages of the pandemic and the community quarantine intervention has shown to have a positive effect in mitigating the spread of the disease. However, when the government attempted to ease the economic burden by implementing the MGCQ in Q3, the R d t increased indicating a faster disease spread in the city. The reimplementation of GCQ in Q4 has reduced the R d t to a value of slightly less than one.
It is noteworthy however that at Q3, a difference between the R d t and R s t can be seen. As Q3 started, the initial conditions of the transmission dynamics were during a spike of cases yielding twice as much β value as Q2, hence as discussed in the results of elasticity indices, the spike in β also doubled the R d t . This implies that during the first data points in Q3, R d t captured the spike and the extreme values of R s t . Nevertheless, since R s t is data and time-dependent, its values decreased several days after day one of Q3 due to the perceived lag effect of the intervention implemented in that period. From here, it is conjectured that discrepancies between R d t and R s t in Q3 can be due to the assumed initial conditions (i.e., condition at day 1 of Q3) of the model. Nonetheless, the discrepancies between R d t and R s t are deemed tolerable.

Conclusion
This paper presented a community-level COVID-19 transmission model incorporating the country's health department quarantine and isolation protocols with Davao City, Philippines, as a case study. Key epidemiological parameters were represented and were analyzed in conjunction with the different levels and quarantine classifications (e.g., ECQ, GCQ, and MGCQ). Qualitative analyses indicate that the formulated mathematical model is well-posed, bounded, and locally asymptotically stable when it is disease-free, and when the disease is endemic. Using nonlinear least squares techniques, the parameters β and r were estimated which were used along with other data-estimated parameters to compute epidemiological measures such as the basic reproduction number, and time-varying reproduction number. Uncertainty and sensitivity analyses were also performed via the LHS/PRCC method.
The result on the reproduction numbers (R 0 and R d t ) and elasticity indices as discussed in the Qualitative analyses subsection implied that the number of secondary infections was driven by the number of infectious populations and how often they interact with the susceptible population. The transmission rate in Davao City was estimated to be highest during ECQ, i.e., at the start of the pandemic (Table 3). This indicated that Davao City may had limited and insufficient resources against the emergence of a pandemic. Hence, there is a need to better equip the healthcare system and the local government units (LGU) for better disaster management and prevention of future pandemics. A multi-sectoral framework, policies, and protocols of how the different agencies could collaboratively work should already be put in place even before the occurrence of any health emergencies. The National, Regional, Provincial, City, and Municipal Disaster Risk Reduction and Management Offices, the Department of Health (DOH), the Department of Interior and Local Government, the LGU, Provincial, City, Municipal, and Barangay Health Offices are some of the core agencies necessary for effective pandemic management. Guides, manuals, and protocols should be collaboratively drafted by these agencies and periodic capacity training for its personnel should be incorporated into its mandate for better preparedness and mitigation. Based on the model results, policies toward the reduction of the transmission rate such as the implementation of containment measures, availability and wearing of protective equipment, and sanitation are some of the key factors to consider.
On the other hand, as discussed in the above subsection on Epidemiological parameters and corroborated by Table 3 concerning the transmission rates, an abrupt change from a strict containment measure to a lax one could negatively affect the control measures implemented especially because of a perceived lag effect. A protocol for a gradual easement of containment measures should be implemented instead. Moreover, since we know that the results with respect to the positivity proportion r corroborated the need to put forth better monitoring and identification protocols which include, but are not limited to, testing and contact tracing capabilities, and travel and mobility surveillance. Furthermore, the results for the severity proportion q suggested the need to allocate funding and to study the appropriate ratio of the hospital carrying capacity to the city population because despite the small number of severely symptomatic COVID-19 patients, the hospitals in Davao City still were overwhelmed. More than allocating funds for equipment and buildings, the production of and adequate remuneration for an appropriate number of nursing and healthcare professionals should be a necessity [56]. The need for this action was further supported by the results with respect to the diseaseinduced death rate μ H as it only critically increased in Q4 when many of the hospitals reached its maximum capacity, and many of the healthcare workers were overwhelmed by the number of hospital admissions. Despite having available hospital beds, the number of attending nurses and doctors was inadequate, and the available ones are already highly exhausted. Government interventions such as the provisions of scholarships, opening of job positions, and offering of competitive remuneration among others should have more investments and be made more accessible to all aspiring healthcare workers.
Overall, the outputs generated from this study are potentially beneficial as a basis in making decisions crucial for impeding the spread of COVID-19 and may be a baseline basis for the protocols against future pandemics. These outputs also provided a quantitative measure of the respective effect of the various measures implemented during the different quarantine classifications. This model could be used for other cities and regions in the Philippines to assess the effects of their respective efforts to combat the disease, which is contextualized on a community level. Moreover, researchers and practitioners must be aware of the limitations of compartmental-based modeling (SIR, SEIR, etc.). This approach generally assumes that the community under infection is homogeneous: each human host infects and undergoes infection in the same manner. A comparison of our analysis with agent-based models [57] and other modelling techniques accounting heterogeneity can reveal additional factors that can refine epidemic dynamics and projections. Since our model assumed data prior to the vaccine rollout in Davao City, a suitable refinement of our model should incorporate community-wide inoculation dynamics. Other possible extensions account for risk groups [58], particularly age groups [59], and untraced and undetected infectious individuals in the context of testing and contact tracing [54,60].
Supporting information S1 File. Supplemental file. This file contains solutions to the mathematical theories described in this paper. (PDF)